#include "slep.h"
#include "hr_time.h"
#include "Logistic.h"
#include "epph.h"

/** We get b=b*diag(weight); and Weight and returns m1
*/
double normalizeWeights(OPTS* opt,double *weight,double *y,double* b){
	double *pos, *neg, ps,ns,denominator,sWeight[2],m1;
	pos = (double*) malloc(sizeof(double) * opt->m);
	neg = (double*) malloc(sizeof(double) * opt->m);
	dcopy(opt->m,y,1,pos,1);
	dcopy(opt->m,y,1,neg,1);
	daxpy(opt->m,1.0,oneVector,1,pos,1);
	ps = ddot(opt->m,oneVector,1,pos,1);
	ps = ps/2;
	ns = ((double)opt->m) - ps;

	if(!hasGot(opt,"sWeight")){
		sWeight[0] = opt->sWeight[0];
		sWeight[1] = opt->sWeight[1];
		denominator = ps * sWeight[0] + ns * sWeight[1];
		daxpy(opt->m,-1.0,oneVector,1,neg,1);

		memset(b,0,sizeof(double)* opt->m);

		dscal(opt->m,0.5*sWeight[0]/denominator,pos,1);
		dscal(opt->m,-0.5*sWeight[1]/denominator,neg,1);

		daxpy(opt->m,sWeight[1]/denominator,pos,1,b,1);
		daxpy(opt->m,-sWeight[0]/denominator,neg,1,b,1);

		memset(weight,0,sizeof(double)* opt->m);
		daxpy(opt->m,1.0,pos,1,weight,1);
		daxpy(opt->m,1.0,neg,1,weight,1);
		m1 = sWeight[0]*ps / denominator;
		dscal(opt->m,1.0-m1,pos,1);
		dscal(opt->m,-m1,neg,1);
		memset(b,0,sizeof(double)* opt->m);
		daxpy(opt->m,1.0,pos,1,b,1);
		daxpy(opt->m,1.0,neg,1,b,1);
	}else{
		memset(weight,0,sizeof(double)*opt->m);
		daxpy(opt->m,-1.0/opt->m,oneVector,1,weight,1);
		denominator = 1.0;
		m1 = ps / opt->m;
	}
	free(pos);
	free(neg);
	return m1;
}


FUNVAL LogisticR(double* A, double* y, double z, OPTS opt){
	if(opt.lFlag==0 && opt.mFlag==0 )
		return LogisticRm1(A,y,z,opt);
	else if( opt.lFlag ==0 && opt.mFlag==1)
		return LogisticRm2(A,y,z,opt);
	else if( opt.lFlag == 1 && opt.mFlag == 1)
		return LogisticRm3(A,y,z,opt);
}




// Matlab function [x, c, funVal, ValueL]=LogisticR(A, y, z, opts)
FUNVAL LogisticRm1(double* A, double* y, double z, OPTS opt){

	// Common for all 3 LogisticRm1, LogisticRm2, LogisticRm3
	int m,n,i,iterStep,bFlag;
	double m1,m2;
	double lambda, lambda_max,rsL2, c, cp, ccp, alphap, alpha, L, beta, sc, fun_s, fun_x, fVal, fValp, l_sum, r_sum, gc;
	//double *weight,*weighty,*b, *ATb,*x, *Ax, *Axp, *xp, *xxp, *s , *As, *aa, *bb, *prob, *funVal, *g, *v, *ValueL;
	double *x,*xp,*xxp,*Ax,*Axp, 
		*b, *weight, *weighty,
		*ATb,
		*s,*As,
		*aa,*bb,
		*prob,
		*g, *v,
		*ValueL,*funVal;
	FUNVAL funret;

	m = opt.m;
	n = opt.n;

	x = (double *) malloc( sizeof(double) * n);
	xp = (double *) malloc( sizeof(double) * n);
	xxp = (double *) malloc( sizeof(double) * n);
	s = (double *) malloc( sizeof(double) * n);
	Ax = (double *) malloc( sizeof(double) * m);
	Axp = (double *) malloc( sizeof(double) * m);
	As = (double *) malloc( sizeof(double) * m);


	g = (double *) malloc( sizeof(double) * n);
	v = (double *) malloc( sizeof(double) * n);

	b = (double *) malloc( sizeof(double) * m);
	ATb = (double *) malloc( sizeof(double) * n);


	weight = (double *) malloc( sizeof(double) * m);
	weighty = (double *) malloc( sizeof(double) * m);


	aa = (double *) malloc( sizeof(double) * m);
	bb = (double *) malloc( sizeof(double) * m);

	prob = (double *) malloc( sizeof(double) * m);

	ValueL = (double *) malloc( sizeof(double) * opt.maxIter);
	funVal = (double *) malloc( sizeof(double) * opt.maxIter);

	// Initialize normalization
	initNormalization(&opt,A);

	m1 = normalizeWeights(&opt,weight,y,b);
	m2 = 1.0 - m1;


	// Starting point
	rsL2 = opt.rsL2;
	if(opt.rFlag==0)
		lambda = z;
	else{
		if(z<0.0 || z>1.0){
			funret.errorCode = 1;
			return funret;
		}
		normalizedmv('T',&opt,A,b,ATb);
		lambda_max = idamax(n,ATb,1);
		lambda = z * lambda_max;
		rsL2 = rsL2 * lambda_max;
	}

	if( opt.init == 2){
		memset(x,0,sizeof(double)*n);
		c = log(m1/m2);
	}
	else{
		if( hasGot(&opt,"x0"))
			dcopy(n,opt.x0,1,x,1);
		else
			memset(x,0,sizeof(double)*n);

		if( hasGot(&opt,"c0"))
			c = opt.c0;
		else
			c = log ( m1/m2);
	}

	normalizedmv('N',&opt,A,x,Ax);

	// The Armijo Goldstein line search scheme
	//------------------------------------------------
	bFlag = 0;
	L= 1.0/m + rsL2;

	// weighty = weight .* y;
	dsbmv('U',m,0,1.0,weight,1,y,1,0.0,weighty,1);

	dcopy(n,x,1,xp,1);
	dcopy(m,Ax,1,Axp,1);
	memset(xxp,0,sizeof(double)*n);
	cp = c;
	ccp = 0.0;

	alphap = 0.0;
	alpha = 1.0;

	for(iterStep = 0; iterStep < opt.maxIter; iterStep++){

		// beta = (alphap-1)/alpha;
		// s = x + beta* xxp;
		// sc = c + beta*ccp;
		beta = (alphap-1.0)/alpha;
		dcopy(n,x,1,s,1);
		daxpy(n,beta,xxp,1,s,1);
		sc = c + beta*ccp;

		// As = Ax + beta ( Ax-Axp);
		dcopy(m,Ax,1,As,1);
		dscal(m,1.0+beta,As,1);
		daxpy(m,-beta,Axp,1,As,1);

		// aa = -y.* (As+ sc);
		dcopy(m,As,1,aa,1);
		daxpy(m,sc,oneVector,1,aa,1);
		dtbmv('U','N','N',m,0,y,1,aa,1);
		dscal(m,-1.0,aa,1);

		// fun_s= weight' * ( log( exp(-bb) +  exp(aa-bb) ) + bb )+...
        //    rsL2/2 * s'*s;
		fun_s = 0;
		for(i=0; i<m;i++){
			bb[i] = maxof2(aa[i],0.0);
			fun_s += weight[i] * ( log(exp(-bb[i])+exp(aa[i]-bb[i])) + bb[i]);
		}
		fun_s += rsL2 * ddot(n,s,1,s,1)/2.0;

		//b = -weighty .* (1-prob);
		for(i=0;i<m;i++){
			prob[i] = 1.0/(1.0+exp(aa[i]));
			b[i] = - weighty[i] * (1.0 - prob[i]);
		}


		// gc =  sum(b); the gradient of c
		gc = ddot(m,b,1,oneVector,1);

		// g = ATb
		normalizedmv('T',&opt,A,b,g);

		// g = g + rsL2*s;
		daxpy(n,rsL2,s,1,g,1);


		// xp = x; Axp = Ax;
		// cp = c;
		dcopy(n,x,1,xp,1);
		dcopy(m,Ax,1,Axp,1);
		cp = c;

		for(;;){
			// v = s - g/L;
			dcopy(n,s,1,v,1);
			daxpy(n,-1.0/L,g,1,v,1);

			// c = sc - gc/L;
			c = sc - gc/L;

			// x = sign(v) .* max(abs(v)-lambda/L , 0);
			for(i=0; i<n;i++)
				x[i] = (v[i]/fabs(v[i])) *maxof2(fabs(v[i])-lambda/L,0.0);

			// v = x-s;
			dcopy(n,x,1,v,1);
			daxpy(n,-1.0,s,1,v,1);

			normalizedmv('N',&opt,A,x,Ax);

			// aa = -y .* (Ax+c);
			dcopy(m,Ax,1,aa,1);
			daxpy(m,c,oneVector,1,aa,1);
			dtbmv('U','N','N',m,0,y,1,aa,1);
			dscal(m,-1.0,aa,1);

			// bb= max(aa,0);
			// fun_s= weight' * ( log( exp(-bb) +  exp(aa-bb) ) + bb )+...
			//    rsL2/2 * s'*s;
			fun_x = 0;
			for(i=0; i<m;i++){
				bb[i] = maxof2(aa[i],0.0);
				fun_x += weight[i] * ( log(exp(-bb[i])+exp(aa[i]-bb[i])) + bb[i]);
			}
			fun_x += rsL2 * ddot(n,x,1,x,1)/2.0;

			r_sum = (ddot(n,v,1,v,1) + (c-sc)*(c-sc))/2;
			l_sum = fun_x - fun_s - ddot(n,v,1,g,1) - (c-sc)*gc;

			if( r_sum <= 1.0e-20){
				bFlag = 1;
				break;
			}

			if( l_sum <= r_sum * L)
				break;
			else
				L = maxof2(2*L, l_sum/r_sum);
		}
		alphap = alpha;
		alpha = ( 1.0 + sqrt(4*alpha*alpha +1.0))/2.0;

		ValueL[iterStep] = L;

		// xxp = x - xp; ccp = c- cp;
		dcopy(n,x,1,xxp,1);
		daxpy(n,-1.0,xp,1,xxp,1);
		ccp = c-cp;

		// funVal(iterStep) = fun_x + lambda * sum(abs(x));
		fVal = funVal[iterStep ] = fun_x + lambda*dasum(n,x,1);
		if(iterStep != 0)
			fValp = funVal[iterStep-1];
		else 
			fValp= funVal[iterStep];

		if(bFlag)
			break;

		if(terminationCondition(&opt,fVal,fValp,x,dnrm2(n,xxp,1),dnrm2(n,xp,1),iterStep))
			break;
	}

	funret.errorCode = 0;
	funret.c = c;
	funret.funVal = funVal;
	funret.totIter = iterStep;
	strcpy(funret.type,"LogisticR");
	funret.x = x;
	funret.ValueL = ValueL;
	return funret;
}




// Reformulated problem + Nemirovski's line search scheme
FUNVAL LogisticRm2(double* A, double* y, double z, OPTS opt){

	// Common for all 3 LogisticRm1, LogisticRm2, LogisticRm3
	int m,n,i,iterStep,bFlag;
	double m1,m2;
	double lambda, lambda_max,rsL2, c, cp, ccp, alphap, alpha, L, beta, sc, fun_s, fun_x, fVal, fValp, l_sum, r_sum, gc, norm_xxp, norm_xp;
	//double *weight,*weighty,*b, *ATb,*x, *Ax, *Axp, *xp, *xxp, *s , *As, *aa, *bb, *prob, *funVal, *g, *v, *ValueL;
	double *x,*xp,*xxp,*Ax,*Axp, 
		*b, *weight, *weighty,
		*ATb,
		*s,*As,
		*s_t,
		*t, *tp,
		*aa,*bb,
		*prob,
		*g, *v, *v_t, *u,
		*ValueL,*funVal;
	FUNVAL funret;

	m = opt.m;
	n = opt.n;

	x = (double * ) malloc( sizeof(double) * n);
	xp = (double * ) malloc( sizeof(double) * n);
	xxp = (double * ) malloc( sizeof(double) * n);
	Axp = (double * ) malloc( sizeof(double) * m);
	Ax = (double * ) malloc( sizeof(double) * m);
	b = (double * ) malloc( sizeof(double) * m);
	ATb = (double * ) malloc( sizeof(double) * n);

	weight = (double * ) malloc( sizeof(double) * m);
	weighty = (double * ) malloc( sizeof(double) * m);

	s = (double * ) malloc( sizeof(double) * n);
	s_t = (double * ) malloc( sizeof(double) * n);
	t = (double * ) malloc( sizeof(double) * n);
	tp = (double * ) malloc( sizeof(double) * n);

	As = (double * ) malloc( sizeof(double) * m);
	aa = (double * ) malloc( sizeof(double) * m);
	bb = (double * ) malloc( sizeof(double) * m);
	prob = (double * ) malloc( sizeof(double) * m);

	g = (double * ) malloc( sizeof(double) * n);
	u = (double * ) malloc( sizeof(double) * n);
	v = (double * ) malloc( sizeof(double) * n);
	v_t = (double * ) malloc( sizeof(double) * n);

	ValueL = (double * ) malloc( sizeof(double) * opt.maxIter);
	funVal = (double * ) malloc( sizeof(double) * opt.maxIter);

	// Initialize normalization
	initNormalization(&opt,A);

	m1 = normalizeWeights(&opt,weight,y,b);
	m2 = 1.0 - m1;


	// Starting point
	rsL2 = opt.rsL2;
	if(opt.rFlag==0)
		lambda = z;
	else{
		if(z<0.0 || z>1.0){
			funret.errorCode= 0;
			return funret;
		}
		normalizedmv('T',&opt,A,b,ATb);
		lambda_max = idamax(n,ATb,1);
		lambda = z * lambda_max;
		rsL2 = rsL2 * lambda_max;
	}

	if( opt.init == 2){
		memset(x,0,sizeof(double)*n);
		c = log(m1/m2);
	}
	else{
		if( hasGot(&opt,"x0"))
			dcopy(n,opt.x0,1,x,1);
		else
			memset(x,0,sizeof(double)*n);

		if( hasGot(&opt,"c0"))
			c = opt.c0;
		else
			c = log ( m1/m2);
	}

	normalizedmv('N',&opt,A,x,Ax);

	bFlag = 0;

	L = rsL2 + 1.0/m ;

	// weighty = weight .* y;
	dsbmv('U',m,0,1.0,weight,1,y,1,0.0,weighty,1);	

	// assigning xp with x, and Axp with Ax
	dcopy(n,x,1,xp,1);
	dcopy(m,Ax,1,Axp,1);
	memset(xxp,0, sizeof(double)*n);
	cp = c;
	ccp = 0;

	alphap = 0.0; 
	alpha = 1.0;

	for ( iterStep = 0; iterStep < opt.maxIter; iterStep++ ){
		// computer search point based on xp and x
		beta = (alphap -1.0)/alpha;
		
		// s= x + beta*xxp; sc = c+ beta*ccp; s_t= t + beta* (t-tp);
		dcopy(n,x,1,s,1);
		daxpy(n,beta,xxp,1,s,1);
		sc = c + beta * ccp;
		dcopy(n,t,1,s_t,1);
		dscal(n,1.0+beta,s_t,1);
		daxpy(n,-beta,tp,1,s_t,1);

		// As = Ax + beta* ( Ax-Axp);
		dcopy(m,Ax,1,As,1);
		dscal(m,1.0+beta,As,1);
		daxpy(m,-beta,Axp,1,As,1);

		// aa = - diag(y) * (As + sc)
		dcopy(m,As,1,aa,1);
		daxpy(m,sc,oneVector,1,aa,1);
		dtbmv('U','N','N',m,0,y,1,aa,1);
		dscal(m,-1.0,aa,1);

		// fun_s= weight' * ( log( exp(-bb) +  exp(aa-bb) ) + bb )+...
        //    rsL2/2 * s'*s;
		fun_s = 0;
		for(i=0; i<m;i++){
			bb[i] = maxof2(aa[i],0.0);
			fun_s += weight[i] * ( log(exp(-bb[i])+exp(aa[i]-bb[i])) + bb[i]);
		}
		fun_s += rsL2 * ddot(n,s,1,s,1);

		//b = -weighty .* (1-prob);
		for(i=0;i<m;i++){
			prob[i] = 1.0/(1.0+exp(aa[i]));
			b[i] = - weighty[i] * (1.0 - prob[i]);
		}

		// gc =  sum(b);
		gc = ddot(m,b,1,oneVector,1);

		// g= A'*b;
		normalizedmv('T',&opt,A,b,g);

		// g = g + rsL2 * s;
		daxpy(n,rsL2,s,1,g,1);

		// xp = x; Axp=Ax;
		// cp = c; tp =t;

		dcopy(n,x,1,xp,1);
		dcopy(m,Ax,1,Axp,1);
		cp = c;
		dcopy(n,t,1,tp,1);

		for(;;){
			// c =  sc - gc/L;
			c = sc - gc/L;

			// u = s - g/L;
			// v= s_t - lambda/L;

			dcopy(n,s,1,u,1);
			daxpy(n,-1.0/L,g,1,u,1);
			dcopy(n,s_t,1,v,1);
			daxpy(n,-lambda/L,oneVector,1,v,1);

			ep1R(x,t,u,v,n);
			
			//v = x - s;
			dcopy(n,x,1,v,1);
			daxpy(n,-1.0,s,1,v,1);

			// v_t = t - s_t;
			dcopy(n,t,1,v_t,1);
			daxpy(n,1.0,s_t,1,v_t,1);

			normalizedmv('N',&opt,A,x,Ax);

			// aa = - diag(y) * (A * x + c * 1Vector)
			dcopy(m,Ax,1,aa,1);
			daxpy(m,c,oneVector,1,aa,1);
			dtbmv('U','N','N',m,0,y,1,aa,1);
			dscal(m,-1.0,aa,1);

			fun_x= 0.0;
			for(i = 0 ; i <m ; i++){
				bb[i] = maxof2(aa[i],0);
				fun_x += weight[i] * ( log( exp(-bb[i]) + exp(aa[i]-bb[i]) ) + bb[i]);
			}
			fun_x += rsL2 * ddot(n,x,1,x,1);

			r_sum = ( ddot(n,v,1,v,1) + ddot(n,v_t,1,v_t,1) + (c-sc)*(c-sc))/2.0;
			l_sum = fun_x - fun_s - ddot(n,v,1,g,1) - (c-sc)*gc;
			
			if(r_sum <= 1.0e-20){
				bFlag=1;
				break;
			}

			if(l_sum <= r_sum*L)
				break;
			else
				L=maxof2(2*L,l_sum/r_sum);
		}
		alphap = alpha;
		alpha = (1.0 + sqrt(4.0*alpha*alpha + 1.0))/2;

		ValueL[iterStep] = L;

		// xxp = x-xp;
		//ccp = c-cp;
		dcopy(n,x,1,xxp,1);
		daxpy(n,-1.0,xp,1,xxp,1);
		ccp = c-cp;

		fValp = fVal = funVal[iterStep] = fun_x + lambda* ddot(n,t,1,oneVector,1);

		if(iterStep !=0)
			fValp = funVal [ iterStep -1];

		if(bFlag)
			break;
		//     norm_xxp=sqrt(xxp'*xxp + norm(t-tp)^2 + (c-cp)^2);
		norm_xxp = ddot(n,xxp,1,xxp,1) + ddot(n,t,1,t,1) + ddot(n,tp,1,tp,1) - 2 * ddot(n,t,1,tp,1) + (c-cp)*(c-cp);
		norm_xxp = sqrt(norm_xxp);
		//     norm_xp=sqrt(xp'*xp + cp^2 + tp'*tp);
		norm_xp = sqrt(ddot(n,xp,1,xp,1)+ cp*cp+ ddot(n,tp,1,tp,1));
		if(terminationCondition(&opt,fVal,fValp,x,norm_xxp,norm_xp,iterStep))
			break;
	}
	funret.errorCode = 0;
	funret.c = c;
	funret.funVal = funVal;
	funret.totIter = iterStep;
	strcpy(funret.type,"LogisticR");
	funret.x = x;
	funret.ValueL = ValueL;
	return funret;

}







FUNVAL LogisticRm3(double* A, double* y, double z, OPTS opt){
	// Common for all 3 LogisticRm1, LogisticRm2, LogisticRm3
	int m,n,i,iterStep,bFlag;
	double m1,m2;
	double lambda, lambda_max,rsL2, c, cp, ccp, alphap, alpha, L, beta, sc, fun_s, fun_x, fVal, fValp, l_sum, r_sum, gc, norm_xxp, norm_xp, gamma, cnew, tao;
	FUNVAL funret;
	double *weight,*weighty,*b, *ATb,*x, *Ax, *Axp, *xp, *xxp, *s , *As, *aa, *bb, *prob, *funVal, *g, *v, *ValueL, *t, *tp, *s_t, *u, *xnew, *Axnew, *tnew, *v_t;

	m = opt.m;
	n = opt.n;

	b = (double * ) malloc( sizeof(double) * m);
	ATb = (double * ) malloc( sizeof(double) * n);

	x = (double * ) malloc( sizeof(double) * n);
	xp = (double * ) malloc( sizeof(double) * n);
	xxp = (double * ) malloc( sizeof(double) * n);
	Ax = (double * ) malloc( sizeof(double) * m);
	Axp = (double * ) malloc( sizeof(double) * m);

	weight = (double * ) malloc( sizeof(double) * m);
	weighty = (double * ) malloc( sizeof(double) * m);


	t = (double * ) malloc( sizeof(double) * n);
	tp = (double * ) malloc( sizeof(double) * n);
	s = (double * ) malloc( sizeof(double) * n);
	s_t = (double * ) malloc( sizeof(double) * n);

	As = (double * ) malloc( sizeof(double) * m);

	aa = (double * ) malloc( sizeof(double) * m);
	bb = (double * ) malloc( sizeof(double) * m);
	prob = (double * ) malloc( sizeof(double) * m);

	g = (double * ) malloc( sizeof(double) * n);
	u = (double * ) malloc( sizeof(double) * n);
	v = (double * ) malloc( sizeof(double) * n);
	tnew = (double * ) malloc( sizeof(double) * n);
	xnew = (double * ) malloc( sizeof(double) * n);
	v_t = (double * ) malloc( sizeof(double) * n);

	Axnew = (double * ) malloc( sizeof(double) * m);

	ValueL = (double * ) malloc( sizeof(double) * opt.maxIter);
	funVal = (double * ) malloc( sizeof(double) * opt.maxIter);

	// Initialize normalization
	initNormalization(&opt,A);

	m1 = normalizeWeights(&opt,weight,y,b);
	m2 = 1.0 - m1;


	// Starting point
	rsL2 = opt.rsL2;
	if(opt.rFlag==0)
		lambda = z;
	else{
		if(z<0.0 || z>1.0){
			funret.errorCode = 1;
			return funret;
		}
		normalizedmv('T',&opt,A,b,ATb);
		lambda_max = idamax(n,ATb,1);
		lambda = z * lambda_max;
		rsL2 = rsL2 * lambda_max;
	}

	if( opt.init == 2){
		memset(x,0,sizeof(double)*n);
		c = log(m1/m2);
	}
	else{
		if( hasGot(&opt,"x0"))
			dcopy(n,opt.x0,1,x,1);
		else
			memset(x,0,sizeof(double)*n);

		if( hasGot(&opt,"c0"))
			c = opt.c0;
		else
			c = log ( m1/m2);
	}

	normalizedmv('N',&opt,A,x,Ax);

	bFlag =  0;

	L = rsL2 + 1.0/m;

	// weighty = weight .* y;
	dsbmv('U',m,0,1.0,weight,1,y,1,0.0,weighty,1);

	gamma = 1.0;

	// xp = x; Axp = Ax; xxp = zeros; cp = c; ccp= 0;
	dcopy(n,x,1,xp,1);
	dcopy(m,Ax,1,Axp,1);
	memset(xxp,0,sizeof(double)*n);
	cp = c;
	ccp = 0.0;

	for(i = 0; i <n; i++)
		t[i]= fabs(x[i]);
	dcopy(n,t,1,tp,1);

	for( iterStep = 0; iterStep < opt.maxIter; iterStep++){
		for(;;){
			if(iterStep!=0){
				alpha = (-gamma+sqrt(gamma*gamma+ 4.0*L*gamma));
				beta = (gamma - gamma*alphap)/(alphap * gamma + alphap*L * alpha);

				//s = x + beta * xxp; s_t = t + beta * (t-tp); sc = c + beta * ccp;
				dcopy(n,x,1,s,1);
				daxpy(n,beta,xxp,1,s,1);
				dcopy(n,t,1,s_t,1);
				dscal(n,1.0 + beta,s_t,1);
				daxpy(n,-beta,tp,1,s_t,1);

				sc = c + beta * ccp;

				// As = Ax + beta * ( Ax - Axp );
				dcopy(m,Ax,1,As,1);
				dscal(m,1.0 + beta, As,1);
				daxpy(m,-beta,Axp,1,As,1);
			}
			else{
				alpha = (-1.0 + sqrt(5.0))/2.0;
				beta = 0.0;
				// s=x ; s_t = t; sc = c;
				dcopy(n,x,1,s,1);
				dcopy(n,t,1,s_t,1);
				sc = c;

				// As = Ax;
				dcopy(m,Ax,1,As,1);
			}
			// aa = - diag(y) * ( A*s + sc)
			dcopy(m,As,1,aa,1);
			daxpy(m,sc,oneVector,1,aa,1);
			dtbmv('U','N','N',m,0,y,1,aa,1);
			dscal(m,-1.0,aa,1);

			fun_s = 0.0;
			for( i = 0; i<m;i ++){
				bb[i] = maxof2(aa[i],0.0);
				fun_s += weight[i] * (log( exp(-bb[i])+exp(aa[i]-bb[i])) + bb[i] );
			}
			fun_s += rsL2 * ddot(n,s,1,s,1)/2.0;

			//b = -weighty .* (1-prob);
			for(i=0;i<m;i++){
				prob[i] = 1.0/(1.0+exp(aa[i]));
				b[i] = - weighty[i] * (1.0 - prob[i]);
			}

			// gc =  sum(b);
			gc = ddot(m,b,1,oneVector,1);

			// g= A'*b;
			normalizedmv('T',&opt,A,b,g);

			// g = g + rsL2 * s;
			daxpy(n,rsL2,s,1,g,1);

			//cnew = sc - gc/L;
			// u = s - g/L;
			// v = s_t - lambda /L;
			cnew = sc - gc/L;
			dcopy( n,s,1,u,1);
			daxpy(n,-1.0/L,g,1,u,1);
			dcopy(n,s_t,1,v,1);
			daxpy(n,-lambda/L,oneVector,1,v,1);

			ep1R(xnew,tnew,u,v,n);

			// v = xnew - s;
			dcopy(n,xnew,1,v,1);
			daxpy(n,-1.0,s,1,v,1);

			// v_t = tnew - s_t;
			dcopy(n,tnew,1,v_t,1);
			daxpy(n,-1.0,s_t,1,v_t,1);

			normalizedmv('N',&opt,A,xnew,Axnew);

			// aa = -diag(y) * ( A*xnew + cnew)
			dcopy(m,Axnew,1,aa,1);
			daxpy(m,cnew,oneVector,1,aa,1);
			dtbmv('U','N','N',m,0,y,1,aa,1);
			dscal(m,-1.0,aa,1);

			fun_x = 0.0;
			for( i = 0; i < m; i++ ){
				bb[i] = maxof2(aa[i],0.0);
				fun_x += weight[i]* ( log( exp(-bb[i]) + exp(aa[i] - bb[i])) + bb[i]);
			}
			fun_x += ddot(n,xnew,1,xnew,1) * rsL2/2.0;

			r_sum = (ddot(n,v,1,v,1) + ddot(n,v_t,1,v_t,1) + (cnew-sc)*(cnew-sc))/2;
			l_sum = fun_x - fun_s - ddot(n,v,1,g,1) - (cnew-sc)*gc;

			if(r_sum <= 1.0e-20){
				bFlag = 1;
				break;
			}

			if( l_sum <= r_sum * L)
				break;
			else
				L = maxof2( 2*L, l_sum/r_sum);
		}
		gamma = L * alpha * alpha;
		alphap = alpha;

		ValueL [iterStep] = L;

		tao = L * r_sum/l_sum;
		if(tao >= 5.0)
			L = L * 0.8;
		// xp = x; x= xnew; xxp = x - xp;
		// Axp = Ax; Ax = Axnew;
		// tp = t; t= tnew;
		// cp = c; c = cnew; ccp = c - cp;

		dcopy(n,x,1,xp,1);
		dcopy(n,xnew,1,x,1);
		dcopy(n,x,1,xxp,1);
		daxpy(n,-1.0,xp,1,xxp,1);
		dcopy(m,Ax,1,Axp,1);
		dcopy(m,Axnew,1,Ax,1);
		dcopy(n,t,1,tp,1);
		dcopy(n,tnew,1,t,1);
		cp = c;
		c = cnew;
		ccp = c - cp;

		fValp = fVal = funVal[iterStep] = fun_x + lambda* ddot(n,t,1,oneVector,1);
		if(iterStep !=0)
			fValp = funVal[iterStep - 1];

		if(bFlag)
			break;

		norm_xxp =ddot(n,tp,1,tp,1);
		norm_xp= sqrt(ddot(n,xp,1,xp,1) + norm_xxp + cp*cp);
		norm_xxp += ddot(n,t,1,t,1) - 2*ddot(n,t,1,tp,1) + (c-cp)*(c-cp);
		norm_xxp = sqrt(norm_xxp);

		if(terminationCondition(&opt,fVal,fValp,x,norm_xxp,norm_xp,iterStep))
			break;
	}
	funret.errorCode = 0;
	funret.c = c;
	funret.funVal = funVal;
	funret.totIter = iterStep;
	strcpy(funret.type,"LogisticR");
	funret.x = x;
	funret.ValueL = ValueL;
	return funret;
}